function F = objfun103(x,feedback,settings)
% F = objfun103(x,feedback,settings)
% cost function of MPC103 which controls mean(R^2) and var(R^2)
% var(R^2) is the variance of 100 samples
% mean(R^2) is calculate using analytical solution of alpha^2 and beta^2

Tab_R_DepRate   = [0.10,0.15,0.20,0.50,1.0];
Tab_R_nu        = [0.0451e-4,1.3915e-05,0.3583e-4,0.1232e-3,0.3939e-3];
Tab_R_sigma2    = [1.0810e-3,6.2341e-3,2.1983e-2,0.0990,0.3007];

if nargin == 2
    varR2_set   = 200;
    Fact_varR2  = 0.5;
    R2_set      = 50;
    Fact_r2     = 0.5;
    P           = 5;   % Prediction steps
    m           = 20;  % Mode
    dt          = 10;
elseif nargin == 3
%     assert(isstruct(settings),'\nsettings should be a struct');
    varR2_set   = settings.varR2_set;
    Fact_varR2  = settings.Fact_varR2;
    R2_set      = settings.R2_set;
    Fact_r2     = settings.Fact_r2;
    P           = settings.P;
    m           = settings.m;
    dt          = settings.dt;
else
    error('\nThere must be two or three input');
end

% assert(isstruct(feedback),'\nfeedback should be a struct');
rho         = zeros(P+1,1);
h           = zeros(P+1,1);
alpha2      = zeros(P+1,m);
beta2       = zeros(P+1,m);

rho(1)      = feedback.rho;
h(1)        = feedback.h;
alpha2(1,:) = (feedback.alpha).^2;
beta2(1,:)  = (feedback.beta).^2;

F = 0;
for i = 1:P
    if Tab_R_DepRate(1) <x(i) <= Tab_R_DepRate(end)
        nu = interp1(Tab_R_DepRate,Tab_R_nu,x(i),'linear','extrap');
        sigma2 = interp1(Tab_R_DepRate,Tab_R_sigma2,x(i),'linear','extrap');
    elseif x(i) <= Tab_R_DepRate(1)
        nu = Tab_R_nu(1);
        sigma2 = Tab_R_sigma2(1);
    else
        nu = Tab_R_nu(end);
        sigma2 = Tab_R_sigma2(end);
    end
    
    for j = 1:m
        temp = sigma2/(2*nu*j^2);
        temp2 = exp(-2*nu*j^2*dt);
        alpha2(i+1,j) = temp+(alpha2(i,j)-temp)*temp2;
        beta2(i+1,j) = temp+(beta2(i,j)-temp)*temp2;
    end
    r2 = sum(alpha2(i+1,:)+beta2(i+1,:))/(2*pi);
    cost = Fact_r2*((R2_set-r2)/R2_set)^2;    
    F = F+cost;
end

tspan = feedback.t:dt:feedback.t+P*dt;
var_r2 = varR2(feedback.alpha,feedback.beta,x,tspan,0.01);
cost_varR2 = sum(Fact_varR2*((var_r2-varR2_set)/varR2_set).^2);
F = F+cost_varR2;